1 Setup

setwd("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/all-batches")
suppressPackageStartupMessages({
    library(data.table)
    library(DESeq2)
    library(gplots)
    library(here)
    library(hyperSpec)
    library(limma)
    library(LSD)
    library(magrittr)
    library(matrixStats)
    library(parallel)
    library(pander)
    library(plotly)
    library(RColorBrewer)
    library(scatterplot3d)
    library(tidyverse)
    library(tximport)
    library(VennDiagram)
    library(vsn)
})
suppressPackageStartupMessages({
    source("~/Git/UPSCb/UPSCb-common/src/R/featureSelection.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/gopher.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/plot.multidensity.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/volcanoPlot.R")
})
lfc <- 0.5
FDR <- 0.01
pal <- c(brewer.pal(8,"Dark2"),1)
pal2 <- brewer.pal(9,"Paired") #require package RColorBrewer
cols <- rainbow(17)
hpal <- colorRampPalette(c("blue","white","red"))(100)
mar <- par("mar")

2 Raw data

2.1 Loading of sample information

  • Read the sample information
samples <- read.csv("~/Git/UPSCb/projects/arabidopsis-nutrition-TOR/doc/samples3.csv")
  • Remove unnecessary samples
samples %<>% filter(!grepl("P11554_1",SciLifeID)) %>% 
    filter(! SciLifeID %in% c("P13406_101",
                        "P14066_128",
                        "P14066_133",
                        "P13406_102",
                        "P14066_131")) %>%
    mutate(Nutrition,Nutrition=relevel(Nutrition,"NPS")) %>% 
    mutate(AZD,AZD=relevel(AZD,"DMSO"))

samples <- samples[order(samples$Timepoint, samples$Nutrition, samples$AZD),]

samples %<>% mutate(Timepoint,Timepoint=factor(paste0("T",Timepoint)))

2.1.1 Load the data

  • Call the data
filenames <- list.files("../Salmon", 
                    recursive = TRUE, 
                    pattern = "quant.sf",
                    full.names = TRUE)
  • Name the data
names(filenames) <- sub("_S.*","",sapply(strsplit(filenames, "/"), .subset, 3))
  • Match data <=> sample list
filenames <- filenames[match(samples$SciLifeID,names(filenames))]
filenames <-filenames[!is.na(names(filenames))]
samples <- samples[match(names(filenames),samples$SciLifeID),]
  • Annotate the samples
samples$Conditions <- factor(paste(samples$Timepoint,samples$Nutrition,samples$AZD,sep="_"))
samples$Batch <- factor(substr(samples$SciLifeID,1,8))

3 Expression data

Read the expression at the transcript level

tx <- suppressMessages(tximport(files = filenames, 
                                type = "salmon", 
                                txOut = TRUE))

summarise to genes

tx2gene <- data.frame(TXID=rownames(tx$counts),
                      GENEID=sub("\\.[0-9]+","",rownames(tx$counts)))
gx <- summarizeToGene(tx,tx2gene=tx2gene)
## summarizing abundance
## summarizing counts
## summarizing length
kg <- round(gx$counts)

Sanity check

stopifnot(all(colnames(kg) == samples$SciLifeID))

3.1 Raw data export

#dir.create(file.path("analysis_Tom","salmon"),showWarnings=FALSE,recursive=TRUE)
#write.table(kg,file="/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/nutrition-unormalised-gene-expression_data.csv")
#save(kg, samples, file = "/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/counts.rda")

3.2 Preliminary validations

3.2.1 Check for the genes that are never expressed

sel <- rowSums(kg) == 0 
sprintf("%s%% (%s) of %s genes are not expressed",
        round(sum(sel) * 100/ nrow(kg),digits=1),
        sum(sel),
        nrow(kg))
## [1] "16.9% (5452) of 32310 genes are not expressed"

4 Data normalisation

For visualization, the data is submitted to a variance stabilization transformation using DESeq2. The dispersion is estimated independently of the sample type

dds <- DESeqDataSetFromMatrix(
    countData = kg,
    colData = samples,
    design = ~ Conditions)
## converting counts to integer mode
dds <- estimateSizeFactors(dds)

4.1 Perform a Variance Stabilizing Transformation for plotting

vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)

Write out

#write.csv(vst,"/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/library-size-normalized_variance-stabilized_data_nutrition.csv")

5 Multivariate analysis

5.1 PCA

Principal Component Analysis on the normalized data * Establishment of the PCA

conditions1 <- factor(paste(samples$AZD,samples$Nutrition,sep="_"))
pc <- prcomp(t(vsd))
percent <- round(summary(pc)$importance[2,]*100);percent
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 
##   46   21   14    6    3    2    1    1    1    0    0    0    0    0    0    0 
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 PC31 PC32 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
## PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40 PC41 PC42 PC43 PC44 PC45 PC46 PC47 PC48 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
## PC49 PC50 PC51 PC52 PC53 PC54 PC55 PC56 
##    0    0    0    0    0    0    0    0
  • Graphical representation of PC1 x PC2
conds <- droplevels(conditions1)
plot(pc$x[,1],
     pc$x[,2],
     xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
     ylab=paste("Comp. 2 (",percent[2],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("bottomright",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("topleft",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

  • Graphical representation of PC2 x PC3
plot(pc$x[,2],
     pc$x[,3],
     xlab=paste("Comp. 2 (",percent[2],"%)",sep=""),
     ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("topleft",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("bottomleft",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

  • Graphical representation of PC1 x PC3
plot(pc$x[,1],
     pc$x[,3],
     xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
     ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("bottomleft",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("bottomright",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

6 Differential expression based on all conditions compared to the T0

6.1 Preparation of the design

sel <- samples$Timepoint %in% c("T0","T6","T24")
suppressMessages(dds <- DESeqDataSetFromMatrix(
    countData = kg[,sel],
    colData = samples[sel,],
    design = ~ Conditions))

6.2 Differential expression analysis

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing

6.3 Variance Stabilising Transformation

6.3.1 Perform a Variance Stabilizing Transformation for plotting

vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)

6.4 Contrasts to T0

The contrast by default is the first one (not Intercept)

resultsNames(dds)
##  [1] "Intercept"                          "Conditions_T24_NP_AZD_vs_T0_T0_0"  
##  [3] "Conditions_T24_NP_DMSO_vs_T0_T0_0"  "Conditions_T24_NPS_AZD_vs_T0_T0_0" 
##  [5] "Conditions_T24_NPS_DMSO_vs_T0_T0_0" "Conditions_T24_NS_AZD_vs_T0_T0_0"  
##  [7] "Conditions_T24_NS_DMSO_vs_T0_T0_0"  "Conditions_T24_PKS_AZD_vs_T0_T0_0" 
##  [9] "Conditions_T24_PKS_DMSO_vs_T0_T0_0" "Conditions_T6_NP_AZD_vs_T0_T0_0"   
## [11] "Conditions_T6_NP_DMSO_vs_T0_T0_0"   "Conditions_T6_NPS_AZD_vs_T0_T0_0"  
## [13] "Conditions_T6_NPS_DMSO_vs_T0_T0_0"  "Conditions_T6_NS_AZD_vs_T0_T0_0"   
## [15] "Conditions_T6_NS_DMSO_vs_T0_T0_0"   "Conditions_T6_PKS_AZD_vs_T0_T0_0"  
## [17] "Conditions_T6_PKS_DMSO_vs_T0_T0_0"

6.4.1 Different effects at T6

6.4.1.1 T6_NPS_DMSO

res <- results(dds,name = "Conditions_T6_NPS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPSDMSO <- rownames(res[cutoffs,])
T6NPSDMSO_Low <- rownames(res[cutoff2,])
T6NPSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9675 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 5305 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4370 genes that are repressed

6.4.1.2 MA plot

DESeq2::plotMA(res)

6.4.1.3 Volcano plot

volcanoPlot(res)

6.4.1.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.5 T6_NPS_AZD

res <- results(dds,name = "Conditions_T6_NPS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPSAZD <- rownames(res[cutoffs,])
T6NPSAZD_Low <- rownames(res[cutoff2,])
T6NPSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 8744 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4633 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4111 genes that are repressed

6.4.1.6 MA plot

DESeq2::plotMA(res)

6.4.1.7 Volcano plot

volcanoPlot(res)

6.4.1.8 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.9 T6_PKS_DMSO

res <- results(dds,name = "Conditions_T6_PKS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6PKSDMSO <- rownames(res[cutoffs,])
T6PKSDMSO_Low <- rownames(res[cutoff2,])
T6PKSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9264 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 5006 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4258 genes that are repressed

6.4.1.10 MA plot

DESeq2::plotMA(res)

6.4.1.11 Volcano plot

volcanoPlot(res)

6.4.1.12 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.13 T6_PKS_AZD

res <- results(dds,name = "Conditions_T6_PKS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6PKSAZD <- rownames(res[cutoffs,])
T6PKSAZD_Low <- rownames(res[cutoff2,])
T6PKSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 8155 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4336 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 3819 genes that are repressed

6.4.1.14 MA plot

DESeq2::plotMA(res)

6.4.1.15 Volcano plot

volcanoPlot(res)

6.4.1.16 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.17 T6_NS_DMSO

res <- results(dds,name = "Conditions_T6_NS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NSDMSO <- rownames(res[cutoffs,])
T6NSDMSO_Low <- rownames(res[cutoff2,])
T6NSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 6540 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 3519 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 3021 genes that are repressed

6.4.1.18 MA plot

DESeq2::plotMA(res)

6.4.1.19 Volcano plot

volcanoPlot(res)

6.4.1.20 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.21 T6_NS_AZD

res <- results(dds,name = "Conditions_T6_NS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NSAZD <- rownames(res[cutoffs,])
T6NSAZD_Low <- rownames(res[cutoff2,])
T6NSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9212 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4764 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4448 genes that are repressed

6.4.1.22 MA plot

DESeq2::plotMA(res)

6.4.1.23 Volcano plot

volcanoPlot(res)

6.4.1.24 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.25 T6_NP_DMSO

res <- results(dds,name = "Conditions_T6_NP_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPDMSO <- rownames(res[cutoffs,])
T6NPDMSO_Low <- rownames(res[cutoff2,])
T6NPDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9581 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 5091 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4490 genes that are repressed

6.4.1.26 MA plot

DESeq2::plotMA(res)

6.4.1.27 Volcano plot

volcanoPlot(res)

6.4.1.28 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.1.29 T6_NP_AZD

res <- results(dds,name = "Conditions_T6_NP_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T6NPAZD <- rownames(res[cutoffs,])
T6NPAZD_Low <- rownames(res[cutoff2,])
T6NPAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 7782 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4303 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 3479 genes that are repressed

6.4.1.30 MA plot

DESeq2::plotMA(res)

6.4.1.31 Volcano plot

volcanoPlot(res)

6.4.1.32 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2 Different effects at T24

6.4.2.1 T24_NPS_DMSO

res <- results(dds,name = "Conditions_T24_NPS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPSDMSO <- rownames(res[cutoffs,])
T24NPSDMSO_Low <- rownames(res[cutoff2,])
T24NPSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 10560 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 5132 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 5428 genes that are repressed

6.4.2.2 MA plot

DESeq2::plotMA(res)

6.4.2.3 Volcano plot

volcanoPlot(res)

6.4.2.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.5 T24_NPS_AZD

res <- results(dds,name = "Conditions_T24_NPS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPSAZD <- rownames(res[cutoffs,])
T24NPSAZD_Low <- rownames(res[cutoff2,])
T24NPSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 8416 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4125 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4291 genes that are repressed

6.4.2.6 MA plot

DESeq2::plotMA(res)

6.4.2.7 Volcano plot

volcanoPlot(res)

6.4.2.8 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.9 T24_PKS_DMSO

res <- results(dds,name = "Conditions_T24_PKS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24PKSDMSO <- rownames(res[cutoffs,])
T24PKSDMSO_Low <- rownames(res[cutoff2,])
T24PKSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9702 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4757 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4945 genes that are repressed

6.4.2.10 MA plot

DESeq2::plotMA(res)

6.4.2.11 Volcano plot

volcanoPlot(res)

6.4.2.12 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.13 T24_PKS_AZD

res <- results(dds,name = "Conditions_T24_PKS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24PKSAZD <- rownames(res[cutoffs,])
T24PKSAZD_Low <- rownames(res[cutoff2,])
T24PKSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 7296 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 3657 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 3639 genes that are repressed

6.4.2.14 MA plot

DESeq2::plotMA(res)

6.4.2.15 Volcano plot

volcanoPlot(res)

6.4.2.16 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.17 T24_NS_DMSO

res <- results(dds,name = "Conditions_T24_NS_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NSDMSO <- rownames(res[cutoffs,])
T24NSDMSO_Low <- rownames(res[cutoff2,])
T24NSDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 4668 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 2244 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 2424 genes that are repressed

6.4.2.18 MA plot

DESeq2::plotMA(res)

6.4.2.19 Volcano plot

volcanoPlot(res)

6.4.2.20 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.21 T24_NS_AZD

res <- results(dds,name = "Conditions_T24_NS_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NSAZD <- rownames(res[cutoffs,])
T24NSAZD_Low <- rownames(res[cutoff2,])
T24NSAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 7202 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 3289 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 3913 genes that are repressed

6.4.2.22 MA plot

DESeq2::plotMA(res)

6.4.2.23 Volcano plot

volcanoPlot(res)

6.4.2.24 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.25 T24_NP_DMSO

res <- results(dds,name = "Conditions_T24_NP_DMSO_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPDMSO <- rownames(res[cutoffs,])
T24NPDMSO_Low <- rownames(res[cutoff2,])
T24NPDMSO_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 9534 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4689 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4845 genes that are repressed

6.4.2.26 MA plot

DESeq2::plotMA(res)

6.4.2.27 Volcano plot

volcanoPlot(res)

6.4.2.28 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2.29 T24_NP_AZD

res <- results(dds,name = "Conditions_T24_NP_AZD_vs_T0_T0_0")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
T24NPAZD <- rownames(res[cutoffs,])
T24NPAZD_Low <- rownames(res[cutoff2,])
T24NPAZD_High <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 8920 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 4629 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 4291 genes that are repressed

6.4.2.30 MA plot

DESeq2::plotMA(res)

6.4.2.31 Volcano plot

volcanoPlot(res)

6.4.2.32 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

7 Comparisons of GOI lists (Venn diagrams)

7.1 GOI of nutrition at T6

  • All GOI
grid.draw(venn.diagram(list(T6NPSDMSO, T6PKSDMSO, T6NSDMSO, T6NPDMSO),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

  • Downregulated genes
grid.draw(venn.diagram(list(T6NPSDMSO_Low, T6PKSDMSO_Low, T6NSDMSO_Low, T6NPDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

  • Upregulated genes
grid.draw(venn.diagram(list(T6NPSDMSO_High, T6PKSDMSO_High, T6NSDMSO_High, T6NPDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

7.2 GOI of AZD at T6

  • All GOI
grid.draw(venn.diagram(list(T6NPSAZD, T6PKSAZD, T6NSAZD, T6NPAZD),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

  • Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6PKSAZD_Low, T6NSAZD_Low, T6NPAZD_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

  • Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6PKSAZD_High, T6NSAZD_High, T6NPAZD_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("Full","No N","No P", "No C")))

7.3 Nutrition and AZD interaction at T6

7.3.1 For Carbon starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T6NPAZD_Low, T6NPDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))

  • Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T6NPAZD_High, T6NPDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))

7.3.2 For Phosphorus starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T6NSAZD_Low, T6NSDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))

  • Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T6NSAZD_High, T6NSDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))

7.3.3 For Nitrogen starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_Low, T6NPSDMSO_Low, T6PKSAZD_Low, T6PKSDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Nit x AZD", "-Nit")))

  • Upregulated genes
grid.draw(venn.diagram(list(T6NPSAZD_High, T6NPSDMSO_High, T6PKSAZD_High, T6PKSDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Nit x AZD", "-Nit")))

7.4 Nutrition and AZD interaction at T24

7.4.1 For Carbon starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_Low, T24NPSDMSO_Low, T24NPAZD_Low, T24NPDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))

  • Upregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_High, T24NPSDMSO_High, T24NPAZD_High, T24NPDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Suc x AZD", "-Suc")))

7.4.2 For Phosphorus starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_Low, T24NPSDMSO_Low, T24NSAZD_Low, T24NSDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))

  • Upregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_High, T24NPSDMSO_High, T24NSAZD_High, T24NSDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Phos x AZD", "-Phos")))

7.4.3 For Nitrogen starvation

  • Downregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_Low, T24NPSDMSO_Low, T24PKSAZD_Low, T24PKSDMSO_Low),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Nit x AZD", "-Nit")))

  • Upregulated genes
grid.draw(venn.diagram(list(T24NPSAZD_High, T24NPSDMSO_High, T24PKSAZD_High, T24PKSDMSO_High),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Ctrl","-Nit x AZD", "-Nit")))

8 Export list of DEGs

write.table(T6NPSDMSO_Low, file = "T6NPSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6PKSDMSO_Low, file = "T6PKSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6NSDMSO_Low, file = "T6NSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6NPDMSO_Low, file = "T6NPDMSO_Low.txt", sep = "\t", row.names = FALSE)

write.table(T6NPSDMSO_High, file = "T6NPSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T6PKSDMSO_High, file = "T6PKSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T6NSDMSO_High, file = "T6NSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T6NPDMSO_High, file = "T6NPDMSO_High.txt", sep = "\t", row.names = FALSE)

write.table(T6NPSAZD_Low, file = "T6NPSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6PKSAZD_Low, file = "T6PKSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6NSAZD_Low, file = "T6NSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T6NPAZD_Low, file = "T6NPAZD_Low.txt", sep = "\t", row.names = FALSE)

write.table(T6NPSAZD_High, file = "T6NPSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T6PKSAZD_High, file = "T6PKSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T6NSAZD_High, file = "T6NSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T6NPAZD_High, file = "T6NPAZD_High.txt", sep = "\t", row.names = FALSE)


write.table(T24NPSDMSO_Low, file = "T24NPSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24PKSDMSO_Low, file = "T24PKSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24NSDMSO_Low, file = "T24NSDMSO_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24NPDMSO_Low, file = "T24NPDMSO_Low.txt", sep = "\t", row.names = FALSE)

write.table(T24NPSDMSO_High, file = "T24NPSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T24PKSDMSO_High, file = "T24PKSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T24NSDMSO_High, file = "T24NSDMSO_High.txt", sep = "\t", row.names = FALSE)
write.table(T24NPDMSO_High, file = "T24NPDMSO_High.txt", sep = "\t", row.names = FALSE)

write.table(T24NPSAZD_Low, file = "T24NPSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24PKSAZD_Low, file = "T24PKSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24NSAZD_Low, file = "T24NSAZD_Low.txt", sep = "\t", row.names = FALSE)
write.table(T24NPAZD_Low, file = "T24NPAZD_Low.txt", sep = "\t", row.names = FALSE)

write.table(T24NPSAZD_High, file = "T24NPSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T24PKSAZD_High, file = "T24PKSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T24NSAZD_High, file = "T24NSAZD_High.txt", sep = "\t", row.names = FALSE)
write.table(T24NPAZD_High, file = "T24NPAZD_High.txt", sep = "\t", row.names = FALSE)

9 GO analysis

9.1 T6 DMSO

9.1.1 Induced genes

GO <- gopher(T6NPSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
## Loading required package: jsonlite
## 
## Attaching package: 'jsonlite'
## The following object is masked from 'package:purrr':
## 
##     flatten
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPSDMSO_High <- GO

GO <- gopher(T6NSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NSDMSO_High <- GO

GO <- gopher(T6PKSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6PKSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6PKSDMSO_High <- GO

GO <- gopher(T6NPDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPDMSO_High <- GO

9.1.2 Repressed genes

GO <- gopher(T6NPSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPSDMSO_Low <- GO

GO <- gopher(T6NSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NSDMSO_Low <- GO

GO <- gopher(T6PKSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6PKSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6PKSDMSO_Low <- GO

GO <- gopher(T6NPDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPDMSO_Low <- GO

9.2 T6 AZD

9.2.1 Induced genes

GO <- gopher(T6NPSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPSAZD_High <- GO

GO <- gopher(T6NSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NSAZD_High <- GO

GO <- gopher(T6PKSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6PKSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6PKSAZD_High <- GO

GO <- gopher(T6NPAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPAZD_High <- GO

9.2.2 Repressed genes

GO <- gopher(T6NPSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPSAZD_Low <- GO

GO <- gopher(T6NSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NSAZD_Low <- GO

GO <- gopher(T6PKSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6PKSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6PKSAZD_Low <- GO

GO <- gopher(T6NPAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T6NPAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T6NPAZD_Low <- GO

9.3 T24 DMSO

9.3.1 Induced genes

GO <- gopher(T24NPSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPSDMSO_High <- GO

GO <- gopher(T24NSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NSDMSO_High <- GO

GO <- gopher(T24PKSDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24PKSDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24PKSDMSO_High <- GO

GO <- gopher(T24NPDMSO_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPDMSO_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPDMSO_High <- GO

9.3.2 Repressed genes

GO <- gopher(T24NPSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPSDMSO_Low <- GO

GO <- gopher(T24NSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NSDMSO_Low <- GO

GO <- gopher(T24PKSDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24PKSDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24PKSDMSO_Low <- GO

GO <- gopher(T24NPDMSO_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPDMSO_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPDMSO_Low <- GO

9.4 T24 AZD

9.4.1 Induced genes

GO <- gopher(T24NPSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPSAZD_High <- GO

GO <- gopher(T24NSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NSAZD_High <- GO

GO <- gopher(T24PKSAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24PKSAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24PKSAZD_High <- GO

GO <- gopher(T24NPAZD_High,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPAZD_High.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPAZD_High <- GO

9.4.2 Repressed genes

GO <- gopher(T24NPSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPSAZD_Low <- GO

GO <- gopher(T24NSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NSAZD_Low <- GO

GO <- gopher(T24PKSAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24PKSAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24PKSAZD_Low <- GO

GO <- gopher(T24NPAZD_Low,background=rownames(vsd)[rowSums(vsd)>0],url="athaliana")
GO$go$name
## NULL
GO$kegg$id
## NULL
GO$pfam$name
## NULL
write.table(data.frame(GO$go$id, GO$go$padj), file = "GO_T24NPAZD_Low.txt", sep = "\t",
            row.names = FALSE)
GO_T24NPAZD_Low <- GO

10 Session Info

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] jsonlite_1.6                vsn_3.54.0                 
##  [3] VennDiagram_1.6.20          futile.logger_1.4.3        
##  [5] tximport_1.14.0             forcats_0.4.0              
##  [7] stringr_1.4.0               dplyr_0.8.3                
##  [9] purrr_0.3.3                 readr_1.3.1                
## [11] tidyr_1.0.0                 tibble_2.1.3               
## [13] tidyverse_1.3.0             scatterplot3d_0.3-41       
## [15] RColorBrewer_1.1-2          plotly_4.9.1               
## [17] pander_0.6.3                magrittr_1.5               
## [19] LSD_4.0-0                   limma_3.42.0               
## [21] hyperSpec_0.99-20180627     ggplot2_3.2.1              
## [23] lattice_0.20-38             here_0.1                   
## [25] gplots_3.0.1.1              DESeq2_1.26.0              
## [27] SummarizedExperiment_1.16.0 DelayedArray_0.12.0        
## [29] BiocParallel_1.20.0         matrixStats_0.55.0         
## [31] Biobase_2.46.0              GenomicRanges_1.38.0       
## [33] GenomeInfoDb_1.22.0         IRanges_2.20.1             
## [35] S4Vectors_0.24.0            BiocGenerics_0.32.0        
## [37] data.table_1.12.6          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       rprojroot_1.3-2        htmlTable_1.13.2      
##  [4] XVector_0.26.0         fs_1.3.1               base64enc_0.1-3       
##  [7] rstudioapi_0.10        affyio_1.56.0          bit64_0.9-7           
## [10] AnnotationDbi_1.48.0   lubridate_1.7.4        xml2_1.2.2            
## [13] splines_3.6.1          geneplotter_1.64.0     knitr_1.26            
## [16] zeallot_0.1.0          Formula_1.2-3          broom_0.5.2           
## [19] annotate_1.64.0        cluster_2.1.0          dbplyr_1.4.2          
## [22] BiocManager_1.30.10    compiler_3.6.1         httr_1.4.1            
## [25] backports_1.1.5        assertthat_0.2.1       Matrix_1.2-18         
## [28] lazyeval_0.2.2         cli_1.1.0              formatR_1.7           
## [31] acepack_1.4.1          htmltools_0.4.0        tools_3.6.1           
## [34] affy_1.64.0            gtable_0.3.0           glue_1.3.1            
## [37] GenomeInfoDbData_1.2.2 Rcpp_1.0.3             cellranger_1.1.0      
## [40] vctrs_0.2.0            preprocessCore_1.48.0  gdata_2.18.0          
## [43] nlme_3.1-142           xfun_0.11              rvest_0.3.5           
## [46] testthat_2.3.0         lifecycle_0.1.0        gtools_3.8.1          
## [49] XML_3.98-1.20          zlibbioc_1.32.0        scales_1.1.0          
## [52] hms_0.5.2              lambda.r_1.2.4         curl_4.2              
## [55] yaml_2.2.0             memoise_1.1.0          gridExtra_2.3         
## [58] rpart_4.1-15           latticeExtra_0.6-28    stringi_1.4.3         
## [61] RSQLite_2.1.2          highr_0.8              genefilter_1.68.0     
## [64] checkmate_1.9.4        caTools_1.17.1.2       rlang_0.4.2           
## [67] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [70] htmlwidgets_1.5.1      bit_1.1-14             tidyselect_0.2.5      
## [73] R6_2.4.1               generics_0.0.2         Hmisc_4.3-0           
## [76] DBI_1.0.0              pillar_1.4.2           haven_2.2.0           
## [79] foreign_0.8-72         withr_2.1.2            survival_3.1-7        
## [82] RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.5          
## [85] crayon_1.3.4           futile.options_1.0.1   KernSmooth_2.23-16    
## [88] rmarkdown_1.18         locfit_1.5-9.1         readxl_1.3.1          
## [91] blob_1.2.0             reprex_0.3.0           digest_0.6.23         
## [94] xtable_1.8-4           munsell_0.5.0          viridisLite_0.3.0